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Abstract 

Antisymmetrized molecular dynamics (AMD) with quantum branching pro- 
cesses is reformulated so that it can be applicable to the collisions of heavy 
nuclei such as 197 Au + 197 Au multifragmentation reactions. The quantum 
branching process due to the wave packet diffusion effect is treated as a ran- 
dom term in a Langevin-type equation of motion, whose numerical treatment 
is much easier than the method of the previous papers. Furthermore a new 
approximation formula, called the triple-loop approximation, is introduced in 
order to evaluate the Hamiltonian in the equation of motion with much less 
computation time than the exact calculation. A calculation is performed for 
the 197 Au + 197 Au central collisions at 150 MeV/nucleon. The result shows 
that AMD almost reproduces the copious fragment formation in this reaction. 

PACS numbers: 24.10. Cn, 02.50.Ey, 02.70.Ns, 25.70.Pq 
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I. INTRODUCTION 



Various kinds of microscopic dynamical models have been developed in order to un- 
derstand the various phenomena in heavy ion collisions in the medium energy region. The 
mean field models, such as the time-dependent Hartree-Fock (TDHF) theory and the Vlasov- 
Uehling-Uhlenbeck (VUU) equation Jl|Q|, are good at the precise description of the single 
particle dynamics in the mean field. On the other hand, the advantage of the molecular dy- 
namics models is, generally speaking, that they can describe the many-body correlation 
which is essential in the fragment formation. 

Nuclear multifragmentation has been a hot topics in these years. It can be regarded as a 
manifestation of the liquid-gas phase transition in the nuclear matter, and we expect that the 
precious information of the finite-temperature nuclear matter in high and low density can be 
obtained by studying the multifragmentation reactions. Although the multifragmentation 
should be related to the property of nuclear matter, the ideally equilibrated nuclear matter is 
not formed in the real reactions, and therefore the studies by using the microscopic dynamical 
models are indispensable. Furthermore the multifragmentation is a good touchstone for the 
microscopic models because it includes the nontrivial mechanism for the fragment formation, 
namely the appearance of new cluster correlations with dynamical symmetry breaking from 
an almost uniform excited matter. 

From the viewpoint of the time- dependent quantum theories which solve the time evo- 
lution of the system from the given initial state to the final state, the multifragmentation 
is not easy to treat, because the final state should be a superposition of a huge number of 
the channel wave functions. There are a huge number of possible ways to decompose the 
total system into fragments. The initial state of the reaction and the individual channels of 
the intermediate and final states may be well described by using rather simple wave func- 
tions, but the total wave function of the intermediate or final state is, of course, too much 
complicated to handle [See Fig. However, since the interference among channels is not 
so important in usual cases, the quantum branchings from a single channel to the super- 
position of many channels can be treated as stochastic branching processes without taking 
account of the interference among channels. Namely, in a practical time-dependent model 
where each channel is described by a rather simple wave function, the time evolution of the 
system should be determined by the successive stochastic quantum branching processes in 
addition to the deterministic time evolution within each channel. The necessary quantum 
branching process varies according to the model because it should depend on how the chan- 
nel wave function is restricted to the simple one. The physical observables are calculated as 
the ensemble average values of the expectation values all over the channels. 

As is well known, TDHF is not suitable for the reactions with many channels such as 
multifragmentation because it does not take account of the quantum branchings mentioned 
above. Although a single Slater determinant may be sufficient for the initial two nuclei and 
the fragments in each final channel, it is far from sufficient for the superposition of the final 
channel wave functions. In such cases, what one can expect by solving the deterministic time 
evolution is, at best, that one of the possible channels appears as the final state. In bad cases, 
however, the obtained final state looks like none of the final channels. The latter may be the 
case for the TDHF application to the multifragmentation of an expanding system because 
the mean field gets weaker and weaker as the system expands widely and therefore the 
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FIG. 1. A schematic picture of the quantum branching processes for multichannel reactions. 

diffusing single particle wave functions will never relocalize to form clusters. It is dangerous 
to apply the mean field concept to dilute system where the system may be branching into 
channels whose mean fields should be different from channel to channel. The VUU equation 
can be regarded as one of the extensions of TDHF, where the two-body collision effect is 
introduced as a term in the equation of the one-body phase space distribution. However, 
since this collision term is a deterministic term to take account of the averaged effect of the 
two-nucleon collisions, the quantum branching effect is not included in the VUU equation. 
The Boltzmann Langevin approach on the other hand, introduced the fluctuation as 

a random term associated with the two-nucleon collisions, and it may be applicable to the 
multichannel reactions such as multifragmentation. However, it is difficult to understand 
that the real origin of the cluster formation were the fluctuation due to the two-nucleon 
collisions, because the two-nucleon collisions are rare in the expanding nuclear matter, while 
the cluster formation should take place even in the ideal situation where the nuclear matter 
is uniformly expanding without initial fluctuation. 

On the contrary to the mean field models mentioned above, the molecular dynamics 
models restrict the channel wave function to an (antisymmetrized) product of wave packets. 
The shape of the wave packets are usually kept fixed and the many body wave function 
is parametrized only by the centroids of the wave packets. The benefit of this restriction 
is that we can avoid encountering the situation where the single particle wave functions 
have been expanded and then the mean field concept does not work any longer. In other 
words, the channel wave function in the molecular dynamics models are restricted so that 
it cannot be a mixture of many channels which should be treated independently rather 
than as a whole in a single wave function. In the quantum molecular dynamics (QMD) 
and the antisymmetrized molecular dynamics (AMD) [^,0, the centroid motion in each 
channel is determined by the equation of motion derived from the time-dependent variational 
principle. In addition to it, the effect of the two-body collisions is introduced as a stochastic 
branching process, which brings the system from a single initial state to many possible final 
channels randomly. Because of these reasons, QMD and AMD are suitable frameworks for 
the fragmentation phenomena. 

Many people have ever tried to extend the molecular dynamics models by generalizing 
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the wave packets p|-[Tl|,|6|, usually by treating the width parameters of the wave packets 
as dynamical variables as well as the wave packet centroids. Although this extension can 
be an improvement for some phenomena UTTf , it is rather a change for the worse in the 
context of multifragmentation. Such extension draws the molecular dynamics models close 
to a mean field model which fails into the pathological situation where many independent 
channels are mixed in a single Slater determinant. In fact, Kiderlen et al. f9| and Chomaz 
et al. JTU] reported that the diffused wave packets never shrink again to form clusters in a 
hot expanding nuclear system. 

In Ref. [12], we took a different way to extend AMD by taking account of the precise 



one-body dynamics without losing the benefit of the molecular dynamics models that the 
channel wave function cannot fall into a mixture of many channels. It was achieved not by 
generalizing the channel wave function but by introducing the wave packet diffusion effect 
as a new quantum branching process. This extended AMD is called AMD-V, since the 
wave packet diffusion effect is calculated with the Vlasov equation jTJj]. When a expanding 
system is calculated by AMD-V, one can imagine that not only the centroids expand but 
also the successive quantum branchings take place due to the wave packet diffusion effect, 
and that the multifragment channels appear stochastically. In fact, we showed in Refs. 
]T2j , [TB|] that AMD-V works very well for the multifragmentation in 40 Ca + 40 Ca reaction at 
35 MeV/nucleon, though the expansion in this case is not spherical but mainly in the beam 
direction like the neck fragmentation. The wave packet diffusion process is also related to the 
nucleon emission rate and the energy carried out by emitted nucleons, which was essential 
for the correct prediction of the excitation energies of the produced fragments. No other 
microscopic models have ever reproduced these fragmentation data so nicely. 

Ohnishi and Randrup take yet another approach to improve the molecular dynamics 
models ||15[|. Based on the idea that the essential part of the multifragmentation is governed 



by the statistical effect, they introduced a fluctuation-dissipation term to the equation of 
motion by hand so as to ensure the correct equilibrium property. Although the good statistics 
is the advantage of their model, there is no microscopic and dynamical background for the 
added fluctuation-dissipation term. We would like to emphasize here that their approach is 
not the unique way to get the quantum statistics in the molecular dynamics. Even though 
we start with the microscopic dynamical consideration, it is possible to get the quantum 
statistics as shown in Refs. fll6| , |PT|l . 

In spite of the fact that interesting high-quality multifragmentation data were published 
for heavy system such as 197 Au + 197 Au collisions [[H]] , no satisfactory explanation by micro- 
scopic dynamical models has been given. This difficulty is due to the essentially quantum 
mechanical feature of multifragmentation. Although AMD-V is one of few realistic models 
that have possibility to reproduce the data, it was impossible so far to apply AMD to heavy 
systems because of the reason of CPU time. The main purpose of this paper is, therefore, to 
give a framework of AMD-V whose numerical calculation is feasible even for 197 Au + 197 Au 
collisions, by introducing an improvement and an approximation to the original AMD-V 
framework. 

The necessary CPU time of the original AMD calculation is proportional to the forth 
power of the mass number of the system. This means that the required CPU time for 
197 Au + 197 Au reactions is about 600 times as much as for 40 Ca + 40 Ca reactions. In order 
to overcome this problem, we introduce in this paper a new approximation for the AMD 
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Hamiltonian which can be evaluated with the CPU time proportional to the third power of 
the mass number. This approximation is called the triple-loop approximation. 

In the original AMD-V calculation of Ref. [T^|, the most time-consuming part was the 



procedure to ensure the energy conservation after the quantum branching process due to the 
wave packet diffusion effect. It was necessary to solve a kind of frictional cooling equation at 
least for several time steps to search the energy conserving point. This procedure becomes 
unnecessary and the framework becomes more transparent in this paper when the wave 
packet diffusion effect is reformulated as a random term in a Langevin-type equation of 
motion, the formal structure of which is similar to the equation by Ohnishi and Randrup 
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This paper is organized as follows. In Sec. II, the framework of the improved AMD-V is 
given. Especially the wave packet diffusion process is formulated as a random term in the 
equation of motion. In Sec. Ill, the triple-loop approximation for the AMD Hamiltonian is 
formulated and some tests for this approximation are given. A demonstrative calculation 
for 197 Au + 197 Au collisions at 150 MeV/nucleon is given in Sec. IV in order to show that 
the AMD-V calculation for heavy system is really possible and it is likely to reproduce the 
multifragmentation data. Section V is devoted for the summary. 



II. FRAMEWORK OF AMD WITH QUANTUM BRANCHING PROCESSES 

In constructing a time-dependent quantum model for medium energy heavy ion collisions, 
one should keep in mind the fact that the initial state branches into a huge number of reaction 
channels in the intermediate states and the final state. It is too difficult to follow the time 
evolution of the total many-body wave function in which the many-body correlations are not 
negligible. Therefore, we treat separately the branching into channels and the time evolution 
within each channel. Approximations such as the mean field theory may be valid within 
each channel, while the interference among the branched channels may be unimportant. The 
independence of the time evolution of each channel should be respected. 



A. Channel wave function and equation of motion 

We describe each channel wave function by an AMD wave function which is a single 
Slater determinant of Gaussian wave packets ||, 



$(Z) = det 



exp{-i/( ri - ^=f + ^} Xa% (j)}, (1) 



where the complex variables Z = {Zj; i = 1, . . . , A} = {Z i(T ; i = 1, . . . , A, a = x,y, z} 
represent the centroids of the wave packets. We take the width parameter v = 0.16 fm -2 
and the spin isospin states Xa t = P T , P I, n T , or n |. 

The AMD wave function ([!]) seems to be very simple, but it is sufficient for the de- 
scription of the ground states of nuclei. For example, the binding energies obtained by the 
frictional cooling method |19|] coincide with the experimental data within the precision of 1 
MeV/nucleon even though the common values of v and Tq (mentioned later) are used for 
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all nuclei pO| . Therefore the initial state of the reaction and the individual channel wave 
functions in the intermediate and final states are well described by the AMD wave functions. 

The time evolution of the wave packet centroids Z within the same channel is determined 
by the time dependent variational principle, 

J (*(Z)|*(Z)) ' {> 

from which one can derive the equation of motion for Z, 

dZ jT OH 

Ciajr with a, t = x, y, z is a hermitian matrix defined by 

Curjr = dZ f dZ 1 °&W Z )W Z ))> W 

and TC is the expectation value of the Hamiltonian after the subtraction of the spurious 
kinetic energy of the zero-point oscillation of the center-of-masses of fragments ||, 

n{z) = (HZ)MZ)) - -2M A + To{A - Nf{z)) - (5) 



^E^W^f - . (3) 



The quantum Hamiltonian 



^ = E^ + E«« ( 6 ) 

j=l ZIV1 i<j 



includes an effective two-body interaction such as the Gogny force PT] which can be density 



dependent, since the model wave function is limited to a single Slater determinant. The 
spurious kinetic energies of the zero-point oscillation of the center-of-mass of the isolated 
fragments and nucleons have been subtracted in Eq. (|5[) by introducing a continuous number 
of fragments Np(Z) ||. Without this subtraction, the Q- values for nucleon emissions and 
fragmentations would not be reproduced. The parameter T is 3fr u/2M in principle but 
treated as a free parameter for the adjustment of the binding energies. Once the zero- 
point kinetic energies have been subtracted, the center-of-mass wave function of an isolated 
fragment (or nucleon) should be regarded as a plane wave. This method, however, only takes 
account of the expectation value of the kinetic energy and ignores its quantum fluctuation. 
This point will be reconsidered in the next subsection. 

For the later formulation, it is convenient to introduce the Poisson bracket {J-",Q} and 
the inner product of the canonical gradients (J 7 , Q), 

IT CI = 1 V ^Lr~ l 9Q 
1,1 xhtt dZ<„ ia ' jT dZt 



_1 V ^-C- 1 — (7) 
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1 ' h^dZ ia i,7 ' jT dZ* 



1 ! r , 

+ h£- r dz irT ^dz; T - [ > 

Then the equation of motion (|3]) can be rewritten as 

Z = {Z,H}. (9) 
On the other hand, the frictional cooling equation can be written as 

Z = \{Z,H} + fi(Z,H), (10) 
for which one can show that the energy expectation value decreases as time 

H = fi(H,H)<0 (11) 

for arbitrary A and /i(< 0). 

The one-body Wigner function for the AMD wave function (|l|) is given by 

/(r,p) =8 £e- 2 C , - z fH«^)B tt B-i ) (12) 

ik 

U = ^ r+ 2^ P ' (13) 

where 

B lk = e z ^5 aiak (14) 

is the overlap matrix of the non-orthogonal single particle wave packets in Eq. (|]). In what 
follows, it is sometimes convenient to introduce a QMD-like approximation, 

/(r,p)^8$>- 2 ' u - w *' 2 , (15) 

i 

by using the physical coordinate W = {W,} defined by 

Wi = £(>/5)„ z i. Qu 'W- ( 16 ) 



The coordinates W can be regarded as physical ones because quantities such as the orbital 
angular momentum and the number of the harmonic-oscillator quanta are written in the 
usual way by using W. Furthermore the coordinates W are canonical coordinates when the 
antisymmetrization among more than two packets is negligible. The physical coordinates 
are very useful in various places of the formulation of AMD-V, while the QMD-like approx- 
imation of Eq. (p~5|) is too poor to be useful in the evaluation of the Hamiltonian as will be 
seen in Sec. III. 
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B. Wave packet diffusion process 



The wave packet shape is not allowed to change in AMD. Therefore the dynamics of 
the single particle wave functions is not so precisely described as in TDHF. However, we 
should not extend AMD to TDHF because TDHF has the pathological problem of the 
spurious coupling of channels. In Ref . |12j , instead of extending the channel wave functions, 
we introduced the precise single particle dynamics into AMD as a new stochastic branching 
process. In this subsection, this process is reformulated as a random term of a Langevin-type 
equation of motion which is more suitable for numerical calculations. 



1. Fluctuation due to the wave packet diffusion 

Each nucleon k in an AMD wave function $(Z(£)) of one of the branches at the time t 
is represented approximately by a Gaussian wave packet in phase space 

6 

f k (x,t) =8exp[-2^(x a -X fca (t)) 2 ], (17) 

a=l 

where we have introduced the 6-dimensional phase space coordinate 

l,...,6} = {v^r,p/2/iv^}. (18) 

The centroid {X ka ; a = 1, ... ,6} stands for the physical coordinate W^. In the usual AMD, 
the time evolution of X ka is derived from the equation of motion while the shape of the wave 
packet is fixed. 

However, more reliable time evolution of the one-body distribution function is given by 



the TDHF equation or the Vlasov equation WA 



dt dp dr dr dp 

Writing the expectation value of the Hamiltonian as TC[f] for a Slater determinant repre- 
sented by /(r, p), one can obtain the single-particle Hamiltonian h by 



h{T,p,t)- 



*/(r,p) 



(20) 

/=/(r,P,t) 



for the AMD wave function <f>(Z(t)) whose Wigner function is given by Eq. (|T|). The time 
derivative of the width and shape of the wave packet 

4at(t) = J t J{ X -~ X Ut)){xb - X kb (t))f k ( X ,t)d 6 X (21) 

can be evaluated based on the Vlasov equation (]TQ|) by using the test particle method or by 
the direct analytical calculation. 

It will be useful to note that the wave packet diffusion a\ ah is mainly determined by 
the curvature of the mean field in h(r, p) in the region around the wave packet k. When 
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the potential is quadratic with the curvature ^muj 2 = 2h 2 v 2 /m, the wave packet diffusion 
effect is exactly zero, which is approximately satisfied for the packets inside the nucleus. On 
the other hand, for the packets near the surface of the nucleus, the potential curvature is 
negative and then the wave packet diffusion effect becomes essential. 

Instead of changing the shape of the wave packet f k , we now give fluctuation 5X ka (t) to 
the centroid X^aif) i n order to introduce the wave packet diffusion effect a\. If we assume 
the white noise for the fluctuation, it should satisfy 



SXUt} = 0, (22) 
5X ka (t)SX kb (t') = [&l] ab (t)S(t - t'). (23) 

The negative eigenvalues of b\ have been replaced by zero because the shrinking of the wave 
packet cannot be treated unless we respect the interference among channels. 

Although the higher moments of the fluctuation can also be calculated with the Vlasov 



equation, we expect that their effect is not important. In our early work []12" |, we took the 
distribution function for the fluctuation £ a = 5X ka , 



P(0 = (1 - c)5(0 + c ^|exp(-2X:ea« a6 a), (24) 



^21 



°* = -Y-m*\ (25) 

where c is chosen to give the correct variance of the fluctuation [Eq. fl23|)l. Since a is of 
the order of 1, a big branching takes place with small probability, while no branching takes 
place in most cases. However, we here take the Gaussian distribution 

p(0 = ^f exp H^^ b )' (26) 

= C '[*l)a b \ (27) 



where d is determined by Eq. (P3|) , because this is more convenient for the numerical calcu- 
lation. In this small fluctuation is given to each centroid at every time step. 

It should be noted that the fluctuation 5X ka {t) is spurious for an isolated wave packet k 
because there is no other packets that can absorb the recoil from the fluctuation. Further- 
more the mean field theory [Eq. (|i~9|) 1 is not necessarily valid for light nuclei with A <; 10. 
We should avoid the situation where the unreliable fluctuation for the packets inside a light 
fragment has a drastic effect on the dynamics such as spuriously breaking the fragment. 
Therefore, by checking the packets in the neighborhood of the packet k, we put 5X ka (t) = 
when 



0(1.75- |Re(Zi-Z fe )|J < 10 (28a) 

i 

and 



^0(l.75-|Re(Z i -Z fe )|)Re(Z i -Z 



< 5. (28b) 
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Although this prescription may make the cooling of light fragments too slow, it is not 
a problem practically because the decay of these fragments can be calculated later by a 
statistical decay code. 

For numerical convenience, we now introduce a small delay time r of the response to the 
fluctuation SXk a - The delayed fluctuation Ek a is obtained by the equation 

-E ka (t) = -5X ka (t) - -E ka (t), (29) 

dt T T 

whose solution is 

3 te (f) = - f 5X ka {t')e- {t - tl)lT dt\ (30) 
r Jo 

by assuming Sfc a (0) = for the initial state. Instead of the original fluctuation 5Xk a {t), this 
delayed fluctuation Ek a {t) is to be added to the centroid Xk a (t). In numerical calculations 
we take r = 5 fm/c, which should be smaller than the important time scales of the reaction. 
Since the fluctuation is smoothened by the averaging over the time r, it can be treated easily 
numerically. 

Some readers may be interested in the difference between the fluctuation introduced by 
Ohnishi and Randrup jn| and that of our present work. In our model, the fluctuations of 
different packets (labeled by k) are independent while the correlations of the phase space 
components (labeled by a and b) of each packet is properly incorporated by Eq. (P3"|). This is 
a natural consequence of the fact that our fluctuation is introduced based on the mean field 
model. On the contrary, Ohnishi and Randrup simply ignores the importance of the phase 
space correlations, while they introduce the correlations among different packets without 
any microscopic or dynamical justification. 



2. Equation of motion and conserved quantities 

The above determined fluctuation Ek a {t) or its complex vector representation Hfc(t) is the 
fluctuation to the physical coordinate W^. In order to put it in the equation of motion, it is 
now necessary to convert it to the fluctuation to the original AMD coordinates Z. For this 
purpose let us introduce a time- dependent one-body hermitian operator dkit) that generates 
the fluctuation Efc(t). The form of 6&(t) is taken as 

.4 



Ok(t) = ■ at)|W J (t))(W J (t)| 



|W,(t))(W,(t)|(y* J (t)-a)}, (31) 



where the stochastic complex parameters {ykj(t); j = 1, . . . , A} are to be determined below, 
and 

(32) 

( / W \ 2-, 

(r|W) ocexp{-z/(r- — =J }. (33) 
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In QMD-like approximation by the use of the physical coordinates, the expectation value of 
this one-body operator is calculated as 



OfcW*) = £< w il 3 *l w i> (34) 

i=l 

= i^(y kj (t) • W* - yfc(f) • WOe-I^W-^l" (35) 



By identifying the physical coordinates W with the canonical coordinates, the stochastic 
parameters {y k j(t)} are determined by the requirement that the one-body operator generates 
the fluctuation H^(t) at the moment t, 



dW* 



(36) 

W=W(t) 



Then the fluctuation for Z should be generated as {Z, O k (t)} by the exact expectation value 
of the same one-body operator, 

° k ^ t >- mzmz)) ' (37) 

A 

d k {t) = Y J o kl {t). (38) 

Before putting the fluctuation {Z, O k (t)} in the equation of motion, we should note the 
fact that the fluctuation violates the conservation lows for the total momentum and the 
total energy. Such conservation lows should be achieved through many-body correlations 
in reality. Since this kind of many-body correlations are beyond the scope of the one-body 
dynamics of the Vlasov equation, it is inevitable to introduce the conservation lows by hand. 
By correcting the fluctuation for the conservation lows, the equation of motion of AMD-V 
is now written as 



Z = {Z, H) + £ lk {Z, O k + J2 a km V m } c 

k=l L m k 

+Hk(Z,'H + Y l Pk m Qm) 1- (39) 

in k - 

The first term in the square bracket is the fluctuation due to corrected for the center- 
of-mass coordinate and momentum conservation, and the second term is the cooling (or 
heating) term to ensure the energy conservation. The parameter 7^ can be regarded as 1 
until its meaning is explained later. 

When the system has been decomposed into several clusters, the fluctuation H fc to a 
packet k in one of the clusters should not affect the packets in the other clusters through 
the conservation lows. In order to ensure this point, we define the cluster which includes 
the packet k, where the clusters are identified by the condition that two packets i and j 
belong to the same cluster if |Zj — Zj| < 1.75. The subscript of the Poisson bracket in 
Eq. (|3TS|) indicates that the centroids of the packets in the other clusters are treated as static 
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parameters. Namely, the packets in the other clusters are omitted in the summation in Eq. 
(0), and C^ 1 is replaced by the inverse matrix of the submatrix of C. In Eq. (|39[) , by using the 
Lagrange multipliers a m , the constraints are introduced for the conserved quantities {P m }, 
which are the three components of the center-of-mass coordinate and the three components 
of the total momentum 

(^E^) = ^E^Vv^, (40) 

i i 

(E^) = E 2 ^ImZ i . (41) 

i i 

Then the Lagrange multipliers should be determined by 

{Vi, O k } Ck + £fP*> V m } Ck a km = 0. (42) 

m 

The method to ensure the energy conservation should be considered carefully, because 
it has more drastic effects than the center-of-mass conservation. The set of the packets 
which can be adjusted in order to cancel the energy violation by is restricted to the 
neighborhood of the packet k defined by 

N fc = {z; |Z; - Z k \ < 2.5 and i G C k and % ^ k}. (43) 

The total energy is restored by the frictional cooling term in Eq. ([39|) with fi k adjusted for 
the conservation. Since this cooling term should not violate the other conservation lows, the 
quantities 

{Q m } = {<E r *>, S»> (En*?*)} (44) 

i i i 

are kept constant by determining the Lagrange multipliers f3 km by 

+ E(^' QmhJkm = o. (45) 
m 

The parameter fi k is then determined by 

^ _ Ok + Em ®km'P m }c k 

(H, W + EmAmSm)N t 

in order to conserve the total energy. 

It should be noted that ji k appear in Eq. ( |39"D only through their summation /i = J2 k 
if the constrains are ignored for simplicity. Since \x is an intensive quantity (which is inde- 
pendent of the size of the system) and it is averaged over many independent fluctuations O k , 
one can replace /i with its averaged value Jx which is a function of the current state Z. Then 



the cooling term of Eq. (39) is formally similar to the dissipation term of the Langevin equa- 



tion that Ohnishi and Randrup proposed to introduce together with the fluctuation term 



15| . However, we use Eq. (J46|) directly without replacing \i k with their averaged values, so 



that the total energy is exactly conserved. Furthermore, in our method, we do not need to 
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evaluate the second derivatives of the Hamiltonian 7i which would be necessary in order to 
directly evaluate the averaged value ft. 

The method of the energy conservation is the most difficult ambiguity of this model 
because it is an effect beyond the mean field theories. The above prescription, therefore, 
intend to achieve the energy conservation with the least modification of the other degrees 
of freedom by moving them in the direction of the canonical gradient of the Hamiltonian. 
However, as discussed in Ref. JT2[], it seems that the adjusted degrees of freedom should 



be restricted to the thermal or single-particle ones in order to avoid the unphysical direct 
energy conversion from the collective energy (such as the incident energy of the heavy ion 
collision) to the single-particle energy of the fluctuation. For this purpose, the monopole 
and the quadrupole moments in the coordinate and momentum spaces 

<E r ^ Err) ( 47 ) 

i i 

are also included in {Q m } when is composed of more than 15 packets. 

For an isolated packet k, we have put 5X ka = 0. However, due to the delay time r, the 
delayed fluctuation H fc may not be zero and should be respected even for an isolated packet. 
Therefore, when N k < 4 with Nk being the number of the element of we search the non- 
isolated wave packet i (N > 4) that is the closest to the packet k, and then Nj U N fc U {i} 
and Ci U {k} are used instead of and C k , respectively, in the above formalism. 

Finally we comment on the necessary correction when the system is near the ground 
state. As already discussed in Ref. ||12|| , the fluctuation is small but not exactly zero even 
for the ground state because of the semiclassical nature of the Vlasov equation and the 
restricted Slater determinant in AMD. Since the fluctuation should be zero in the ground 
state, a reduction factor j k is introduced in Eq. (^) in order to cancel the fluctuation only 
near the ground state. By noting that the cooling term becomes zero for the ground state, 
a measure of the difference from the ground state is introduced by 

D k = 6 (7i, W + EftmQ») (48) 

UJV fc J 'cons m 

where N cons denotes the number of the constrained quantities {Q m }, and therefore 6N k — 
N cons is the number of the free degrees of freedom for the energy adjustment. The reduction 
factor 7^ is then taken as 

Ik = 1 1 (49) 

yfl + (/ifc//i fc) 2 

1200 1 5 fm/c f D k \ 3 

^ 6iV fc -iV cons V r U.lMeV/(fm/c)J ' [M) 

so that the coefficient for the cooling term 7fc|/ifc| does not exceed the upper limit /zo/c- The 
purpose of the dependence of fi 0k on (6N k — iV cons ) and r is to make this reduction effect 
independent of the choice of the neighborhood [Eq. (|43|) 1 and the delay time r. With 
this parameterization, the fluctuation is reduced to zero in the ground state, while there is 
almost no reduction soon after hard two-nucleon collisions in the example of 19 A.u + 19 A.u 
collisions shown in Sec. IV. 
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3. Energy fluctuation of emitted packets 



As already discussed, we have subtracted the zero-point kinetic energies of isolated pack- 
ets in Eq. (|]). Therefore the emitted packet k should be regarded as a plane wave of the 
momentum P k = 2h^ s /ulmZk■ This is convenient than treating it as a Gaussian packet 
with a momentum spread because the nucleons in the final channels are usually observed 
experimentally as momentum and energy eigen states and we should ensure the momentum 
and the energy conservations in each cannel. 

Let us consider the case where a wave packet k would be emitted as a whole with the 
momentum centroid value Pofc when the zero-point kinetic energy were not subtracted from 
the Hamiltonian. The momentum of this nucleon is 

Pk = Pofc + q, (51) 



where q is a random number of the Gaussian distribution with (q) = and (qaQr) = f^ud, 
for a, r — x, y, z. The kinetic energy of this wave packet is then 



(XT 



E = P ° fc I P ° fc ' q I q2 (52) 
k 2M M 2M " 1 ' 

When the expectation value of the third term (q 2 )/2M = 3h 2 u/2M is subtracted from the 
Hamiltonian, the added term in Eq. (|5p acts as a repulsive force to this packet. Then it will 
be emitted with the momentum P k that satisfies 

l Ek ) = H = Zjk + (53) 
v 1 2M 2M 2M K ' 

Namely the momentum P k is larger than the true centroid Pofe while the energy expectation 
value does not change because of the absence of the momentum spread when the zero-point 
energy is subtracted. 

This prescription, however, has a shortcoming that it takes account of only the expec- 
tation value of the kinetic energy and ignores its fluctuation. For a preequilibrium nucleon 
in high energy collisions, Pofc may be so large that the fluctuation of the second term of 
Eq. ( |5*^D may play some role though its expectation value is zero. In order to take account 
of this kind of energy fluctuation, we now introduce a random process when each packet is 
emitted. By neglecting the difference of the direction of Pofc and P&, the right amount of 
the energy fluctuation can be produced by changing the momentum as 

Pfc -> {P k ~ $Pk + xA Pk )P k /P k , (54) 

where x is a random number taken from the normal distribution with (x) = and (x 2 ) = 1, 
and 



Sp k = Pk ~ \JPl - h 2 ", (55) 
Ap k = h^pu (56) 

Not only the fluctuation (Ap k ) is given but also the average value of P k is decreased by Sp k 
so that the energy expectation value does not change by this random process. It can be 
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introduced as a new term in the equation for ( p9| ) which is put at the moment when the 
packet k is isolated (Nk = 0) for the first time. The momentum in the above discussion 
should be understood as the relative momentum between the emitted nucleon and the parent 
nucleus. The total energy conservation is achieved by adjusting other degrees of freedom of 
the parent nucleus just in the same way for the fluctuation due to the wave packet diffusion 
effect. 

It should be emphasized here that the above prescription is taken when the wave packet 
is emitted as a whole with the considerably high momentum such as in the early stage 
of high energy collisions. On the other hand, when a low energy nucleon evaporates from a 
nucleus, what comes out of the nucleus is only a high momentum component of the packet as 
we discussed in Refs. [|T^,|l7],|!H • Therefore, the true momentum spread around should be 



small and the above prescription should not be taken in this case. In order to continuously 
connect these two extremes of high and low energies, Eqs. (|55|) and ( |56|) are replaced by 



Sp k = P k - Vmax(p2 - h 2 u, P fc 2 /4), (57) 



Ap k = pP k 5p-5p 2 } (58) 
and the momentum is stochastically changed by Eq. (|54]) . 



C. Two- nucleon collision process 

The combination of the deterministic equation of motion and the quantum branching 
process due to the wave packet diffusion effect is essentially equivalent to the mean field 
theory, such as TDHF, for the short time evolution of a channel wave function. However, 
in medium and high energy collisions, there should be the effect of the residual interaction 
which brings a Slater determinant to a superposition of many Slater determinants. This 
effect is introduced as the stochastic two-nucleon collision process. 

In most molecular dynamics models the stochastic two-nucleon collision process 

has been introduced as the process to cause such branchings. In AMD ||, two-nucleon 
collisions are introduced by the use of the physical coordinates W defined by Eq. ([Top. 
When the physical positions of two nucleons get close, their physical momenta are changed 
randomly according to the differential cross section in a similar way to QMD ||f|]. The 
energy-dependent collision cross section may be modified due to the medium effect which 
can be taken into account as the density dependence of the cross section. The Pauli blocking 
is automatically introduced because of the existence of the Pauli-forbidden region in the 
physical coordinate space ||. 



III. TRIPLE-LOOP APPROXIMATION OF AMD HAMILTONIAN 

Until recently, the application of AMD and AMD-V was limited to relatively light systems 
with the total mass number A < 100, because the CPU time proportional to A 4 is necessary 
for the evaluation of the interaction term in the AMD Hamiltonian, 
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V = 2 X^WiM W< ~ ( Pi ( Pk)Bki' B i3 1 > ( 59 ) 

where y9j are the single particle wave functions in Eq. (|1|) and = (tpi\(pj). In order 
to apply AMD-V to heavy systems such as 197 Au + 197 Au collisions, we now introduce an 
approximation for the AMD Hamiltonian. 

First of all, by using the one-body Wigner function of Eq. (0), the two-body interaction 
term (E% can be rewritten bilinear form of /, 

V = f-v-f 

_ rdrdpdr'dp' , 

= J (2-Kh) 6 " ' p ^ r ' p; r ' p ^ r ' P ^ ' 

where v includes the direct term and the exchange term, 

#(r,p; r ',P') = \ v ( r ~ r ') 

--<%-!■') J dse- t{p - p ' ys/h v(s). (61) 

The spin and isospin degrees of freedom should be implicitly understood. 

The Wigner function is now approximated by a sum of 3A Gaussian functions, 

3 4 

/( U )^nu) = £c P / P G (u), (62) 
P =i 

/ p G (u) = 8e" 2|u - Wp|2 . (63) 

The centroids of f9 are chosen as 
v 

Wj p = i 

Zi + i(Wi-Z<) p = A + z (64) 
Zi-i(W 4 -Z<) p = 2A + i 
for z = 1, . . . , A and p — 1, . . . , 3^4, 

so that the packets cover the important phase space region efficiently In order to get a good 
approximation, the coefficients {c p } are determined by the condition 

fp-v-f" = fp-v-f forp=l,...,3A (65) 

which means that the mean field v ■ f averaged around the phase space point w p should not 
change when the exact Wigner function / is replaced by the approximated one /'". This 
condition is just a linear equation system for {c p }, 

3A 

J2A pq c q = b p , (66) 

q=l 

with 
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4* = (67) 
K = / P G -v ■ f. (68) 

The approximated value of V is then obtained by 

V^V"^f"-v-f" = Y.c P A pq c q . (69) 

pq 

In rare cases, the matrix A pq becomes close to a singular matrix. Then the absolute 
values of c p are large, and /'" oscillates violently in the phase space so as to satisfy Eq. (|65|). 
Since the intent of Eq. (|65| ) is to reproduce / by /"' in the important phase space region 
by requiring /'" to be identical to / around the phase space points {w p }, the resultant 
/'" should be a smooth function for the consistency. Therefore the above equation ( |B"6] ) is 
slightly modified to 

E(^ 2 + ^ = £AA + ^ (70) 

q q 
where e is a small parameter, and 

n f 1 forp = 1,...,A , . 

p \ forp = A + 1,...,3A ■ { ' > 

This modified equation is the condition to minimize the quantity 

E(/ P G • « • f" -f?-v- ff + e 2 £(c P - elf (72) 

p p 

with respect to the coefficients {c p }, so that we can avoid the situation where c p deviates 
from the normal value c p very much. 

It can be easily seen that the approximated interaction V'" can be evaluated with the 
CPU time proportional to A 3 which is necessary for evaluating {b p } by Eq. (|68| ) and for 
solving Eq. flTTj). We can also use a similar approach to approximate the derivatives of V 
with respect to the coordinates Z. The required CPU time is also proportional to A 3 . 

The above formalism can be applied for the density-dependent zero-range force with a 
little extension. The forces like the Gogny force and the Skyrme force have the density- 
dependent term 

{)(r,p;r / ,p / ) = ^o[p(r)] CT 5(r-r'). (73) 

The coefficients {c p } are determined in the same way as above by using Eqs. (|67j), (|68| ) and 
(ffO|) but by replacing the density p(r) in Eq. ([73|) by a constant p . The result {c p } does 
not depend on the value of p . Then the approximated value of V is obtained by 

V ^ V'" = W^V%Vc g , (74) 

pq V HO / 

where p p is a smoothed density around the point Re w p defined by 
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t [fm/c] t [fm/c] 

FIG. 2. Tests of the triple-loop approx- FIG. 3. The same as Fig. but for 

imation (diamonds) along the slow frictional the tests along the dynamics of 40 Ca + 40 Ca 

cooling path for 12 C and 197 Au nuclei, com- collisions at 35 MeV/nucleon. 

pared with the exact values (solid line). The 
results of the QMD-like approximation (dot- 
ted line) are also shown. The expectation 
value of the two-body part of the Gogny 
force is shown as a function of time. 
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c q e 



-fi\ Rewp— Rewq| 2 



(75) 



The parameter /i = | is chosen so as to give a good approximation. 

This triple-loop approximation is tested under various circumstances. Figure |2| shows the 
test along the slow frictional cooling path [Eq. ( |T0"D with A = 1 and [i = —0.25] for two nuclei 
12 C and 197 Au. The randomly excited initial nuclei at t — fm/c are cooled down to the 
ground states at t ~ 300 fm/c. The exact expectation value V/A of the density-independent 
two-body part of the Gogny force is shown by a solid line, while the approximated value 
V" I A with the triple- loop approximation is shown by a diamond for each t. The dotted line 
shows the result of the QMD-like approximation where the expectation value is evaluated 
by using the approximated Wigner function of Eq. flT5|), 



18 




oo 



FIG. 4. Examples of the time evolution of the density projected onto the reaction plane from 
t = fm/c to t = 120 fm/c for central 197 Au + 197 Au collisions at 150 MeV/nucleon. The size of 
the shown area is 80 fm x 80 fm. 

v" = EEV (76) 

p=l q=l 

Compared to the too bad result of the QMD-like approximation, the triple-loop approxima- 
tion always gives a good result within the error of about 1 MeV/nucleon. Figure |3| shows 
the similar information for two events of 40 Ca + 40 Ca collisions at 35 MeV/nucleon. The 
event of the upper part is a peripheral collision and the event of the lower part is a cen- 
tral collision. The triple-loop approximation again gives a sufficiently good result of the 
error within 1 MeV/nucleon. We have also found that the triple-loop approximation for the 
density-dependent force has the precision similar to the density-independent force. 




IV. AN APPLICATION TO Au + Au COLLISIONS 

The 197 Au + 197 Au collisions are calculated by AMD-V for the incident energy 150 
MeV/nucleon and the impact parameter < b < 1 fm. One of the interesting aspects 
in this reaction is the copious formation of the intermediate mass fragments (IMFs) with 
Z > 3 from a strongly expanding system as observed in the experiment of Ref . |H| . On the 



other hand, only few IMFs are produced in the dynamical QMD calculation [18|. Therefore 
it is an important theoretical problem to find out how the copious fragment formation is 
understood in the dynamical framework. 

In the AMD-V calculation presented here, the Gogny force [21] is adopted as the effective 



interaction. It corresponds to the soft equation of state with the incompressibility K = 228 
MeV and the appropriate momentum dependence of the mean field. The expectation value 
of the Hamiltonian is evaluated by using the triple-loop approximation described in the 
previous section. The ground state of the 197 Au nucleus is obtained by the frictional cooling 
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FIG. 5. Calculated charge distribution (histogram) in central 197 Au + 197 Au collisions at 150 



MeV/nucleon, compared with the experimental data (diamonds) of Ref. [18|. The error bars show 
the estimated statistical error of the calculated results. 



method, and it has the reasonable binding energy E/A = 7.4 MeV and the root mean square 
radius (r 2 ) 1//2 = 5.5 fm, while the experimental data are 7.9 MeV and 5.3 fm, respectively. 
The adopted two nucleon collision cross section and the angular distribution are the same 
as those of Ref. ||20|| . Around the two-nucleon collision energy ~ 150 MeV, which is 
important in the present reaction, the pp and nn cross section is the same as the free cross 
section (25 mb). The pn cross section is the same as the free cross section (40 mb) at 
zero-density, but it is reduced to about 30 mb for p > p as the medium effect. 

The produced fragments in the dynamical AMD-V calculation are generally excited and 
their decay is calculated by a statistical model. At every 15 fm/c in the dynamical AMD-V 
calculation, the fragments are identified by linking the two-nucleon pairs with |Z — Zij\/^/u < 
5 fm. With this condition, the identified fragments are well separated spatially in most cases. 
A fragment with the mass number 5 < A < A cr is thrown to the statistical decay code 
directly if its mass number before 15 fm/c was also 5 < A < A CT . Namely, the statistical 
decay of each primordial fragment is calculated when the waiting time £ wa it has passed since 
its mass first became smaller than A CT . The parameters are t wa i t = 22.5 fm/c on the average, 
and A cr is chosen to be 25. However, the dependence of the results on these parameters 
is found to be small. When we take A CT = 20 or 30, or £ wait = 37.5 fm/c on the average, 
the change of the IMF multiplicity is a few percent at most. The adopted statistical decay 
code j|] is based on the sequential binary decay model by Puhlhofer p2 |, but it also takes 
account of the emission of composite particles not only in their ground states but also in 
their excited states with the excitation energy E* < 40 MeV. 

Figure |4] shows the time evolution of the density projected to the reaction plane for two 
events. The total system is once compressed and then expands rather rapidly. From the 
expanding matter, a lot of IMFs are produced. The multiplicity of the primordial IMFs 
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TABLE I. Multiplicities of various particles in central 197 Au + 197 Au collisions at 150 
MeV/nucleon. 



Experiment p| AMD-V 



neutron 

proton 

deuteron 

triton 

3 He 

4 He 

IMF 



92.6 
26.1 
18.6 
17.2 
5.7 
21.0 
10.4 



120.6 
56.8 
14.7 
8.8 
2.3 
16.3 
8.7 



is about 16, and about half of them are to disappear by the statistical decay. Although 
the stopping seems to be strong and the expansion is almost isotropic, the mixture of the 
projectile and the target is not complete. More wave packets of the projectile origin come 
out to the forward direction than to the backward direction. 

In Fig. m the calculated charge distribution is compared with the experimental data 



The calculated result, shown by the solid histogram, reproduces the data very well at least 
in the logarithmic scale. The multiplicities of various particles are compared to the data 
in Table (D. The large IMF multiplicity of the experimental data Mmf = 10.4 is almost 
reproduced by the calculated value Mimf = 8.7, though it is slightly smaller than the data. 
This underestimation is due to the underestimation of the Be and Li multiplicities. We also 
notice in Table [I] that the calculated multiplicities of light particles with 2 < A < 4 are too 
small and the nucleon multiplicity is too large. 

It is useful to consider the gas and liquid parts separately. Here the gas part is composed 
the freely moving nucleons and light particles which are usually emitted after hard two- 
nucleon collisions, and the liquid part is composed of the IMFs which is still bound well by 
the mean field. In order to properly describe the fragment formation, the separation of the 
mass and the energy to the gas and liquid parts is essential as well as the dynamics of the 
liquid part itself. The good reproduction of the IMF multiplicity and the charge distribution 
for Z > 5 suggests that AMD-V describes these aspects very well. On the other hand, the 
failure in the light particle multiplicities can be regarded as a problem in the dynamics of 
the gas part, namely the coalescence of particles in the gas part. 

The light particle and IMF multiplicities are much better reproduced when AMD-V is 
augmented by the coalescence of nucleons and light particles, as will be shown in another 
paper. In this paper, we just mention why the coalescence is not properly treated in AMD- 
V and should be added to AMD-V as an augmenting process. In the medium and high 
energy collisions like the present reaction, a lot of nucleons are emitted. Even though these 
nucleons have almost no correlations among them after hard two-nucleon collisions, a pair 
of a proton and a neutron can form a deuteron when these two nucleons are accidentally 
close to each other in the phase space. In order to correctly predict the probability of the 
coalescence of uncorrelated nucleons, it is necessary for AMD-V to have the correct phase 
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space volume for the bound deuteron state. The phase space volume for the wave packet 
centroids is important because the dynamics is governed by the apparently classical equation 
of motion for the centroids. Since the deuteron is a loosely bound system with a single 
bound state, the phase space volume in AMD is much smaller than the correct quantum 
phase space (2irh) 3 . The deuteron and nucleon yields should therefore be underestimated 
and overestimated, respectively. A large part of tritons and 3 He may also be produced by the 
coalescence mechanism of three nucleons (or a nucleon and a deuteron), and therefore the 
present calculation naturally underestimates their multiplicities. Furthermore, we should 
note that the intrinsic bound states of Li and Be isotopes have the cluster structure of light 
composite particles such as a, t and 3 He, with the small binding energies between them. The 
bound phase space volume in AMD is likely to be smaller than the correct quantum phase 
space, and it is natural that AMD-V underestimates the coalescence of the light composite 
particles to produce Li and Be isotopes directly. 

V. SUMMARY 

The quantum branching processes are essential in the molecular dynamics models in order 
to properly describe the multichannel reactions such as the multifragmentation in heavy ion 
collisions. In addition to the two-nucleon collision process which has been recognized as 
an important process, AMD-V takes account of the wave packet diffusion as a stochastic 
branching process rather than as the shape change of the single particle wave packet like 
in TDHF or the Vlasov equation. AMD-V and the Vlasov equation are equivalent with 
respect to the infinitesimal time evolution of a single Slater determinant without two-nucleon 
collisions, except that AMD-V ignores the interference among the channels and avoids the 
spurious channel correlations. 

In this paper, we reformulated AMD-V in two points so that it is applicable even to 
heavy systems such as 19 A.u + 19 A.u collisions. Fist, the fluctuation due to the wave packet 
diffusion was formulated as a stochastic term in the equation of motion for the wave packet 
centroids. A small Gaussian fluctuation is given to each packet at every time step, instead of 
a big displacement once in a while in the previous framework. This reformulation decreases 
the numerical labor because it simplifies the energy conservation procedure. Second, a new 
triple-loop approximation was introduced for the expectation value of the Hamiltonian with 
respect to the AMD wave function. With this triple-loop approximation, the expectation 
value can be evaluated with the numerical operations proportional to A 3 instead of A A in 
the exact calculation, where A is the mass number of the total system. The error of this 
approximation is about 1 MeV/nucleon at most, and therefore it is useful for the study of 
heavy ion collisions. 

The reformulated AMD-V was applied to 197 An + 19 A.u central collisions at 150 
MeV/nucleon. We adopted the Gogny force as the effective interaction. The calculation 
reproduces the qualitative feature of the experimental data that a lot of fragments are 
produced from the radially expanding system. The large IMF multiplicity was almost re- 
produced by AMD-V quantitatively. This result therefore suggests that AMD-V works well 
for the aspects related to the fragment formation, such as the large energy carried out by 
light particles, the collective expansion, and the appearance of the cluster correlation in the 
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expanding system. However, we found that the nucleon multiplicity is strongly overesti- 
mated and the other light particle multiplicities are underestimated. This should be due to 
the problem of AMD-V in the description of the coalescence of nucleons and light particles 
which is beyond the scope of the current version of AMD-V or other usual mean field theo- 
ries. We will show in another paper how we can incorporate the coalescence to AMD-V and 
that the coalescence process improves the reproduction of the data, including the reactions 
with higher energy. Furthermore, it is also an interesting subject in progress to study the 
fragment formation in relation to the equation of state of the nuclear matter in high and 
low density, and also with isospin asymmetry. 
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